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Digital methods of high precision have been developed for the calculation of electric fields and tra- 
jectories in electrostatic lenses and a computer program was written to apply these methods to the 
two-tube lens. The increased precision results from the use of nine-point formulas in the relaxation 
calculation of potentials in place of previously used five-point formulas and from the use of an improved 
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1. Introduction 

As part of a general program to determine precise 
properties of electrostatic lenses, highly precise digital 
methods were developed for the calculation of electric 
fields and trajectories in electrostatic lenses and a com- 
puter program was written to apply these methods 
to the two-tube lens. 1 Our methods are easily general- 
izable to a wide class of electrostatic lenses. 

The two-tube lens was chosen to test the digital 
methods since considerable data are available for this 
lens. Comparison of our results with the existing data 
shows that our methods have increased the accuracy 
of the calculation of lens properties and makes it pos- 
sible, for the first time, to determine all of the third- 
order aberration coefficients of the two-tube lens for 
meridional trajectories. 

2. Calculation of the Potentials 

2.1. Statement of the Problem 

Since a two-tube lens is cylindrically symmetric, 
the potential distribution inside the cylinders can be 
obtained by solving Laplace's equation in cylindrical 
coordinates (see fig. 1) 
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* Present address: Istituto di Fisica, Universita di Bari. Bari. Italy. 

1 A detailed description of the computer program will be published separately. 



giving as boundary conditions the potentials on a 
contour starting and ending on the axis. Let V\ and 
V2 be the potentials on the two cylinders, D the 
inside diameter, and S the size of the gap between 
the two cylinders. The following boundary conditions 
were used : (see fig. 1) 



AB: 


V=V l 


BC: 


v=v x 


CD: 


v=v, 


DE: 


V = linear variation from V x to i(Ki + V 2 ) 


EF: 


V=\(V^V- i ) 




F(0,0) 
^s-O.ID 



FIGURE 1. Schematic of the two-tube electrostatic lens, showing 
r, z, coordinate system and boundary used in the relaxation c<d- 
cu lotion of potentials. 



AF = BC = \\.W. CD= 1.25D. 
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In the actual calculations the distance AF is typically 
3.5 D. It is sufficient to calculate the potentials only 
inside ABCDEFA since cylindrical symmetry requires 

V(-r,z) = V{r,z) 

and the symmetry about the midplane of the lens 
requires 

V{r,-z)=V l +V 2 -V{r,z). 

In addition, once the solution V (r, z) is obtained for 
the case V\ = 0< V 2 = 1, then the solution for any other 
potentials V u V 2 on the tubes can be obtained from 
the formula 

V(r,z) = V l +(Vt-V l )Vo(r,z) 

2.2. Method of Solution 

Laplace's equation (eq 1) was solved by the relaxa- 
tion method on a matrix of equally spaced points [l]. 2 
The domain ABCDEFA was divided into a network 
of equal sides h. Let i be the index running along the 
z axis, and j the index along r. Let Vij be the potential 
at mesh point (£, j). Then the first order difference 
form of eq (1) is: 

^ +1 -2F / , j +F / , / _ 1 + F / , 7+1 -F / , ; _ l 



h- 



+ 



2hr 



Vi + uj-2y u +Vt. 



i,j . 



h 1 







(2) 



From eq (2) the relationship between V;j and the 
four adjacent points is easily obtained: 



Vij- 



Vi,j + i + Vi,j-i + Vi+i.j+Vi-i t j 



^(Vij+i-Vij^) (3) 

Here we assume that rj is measured in units of h. 
On the axis where 



IdV 

r dr \dr 2 J r = 



(4) 



eq (3) becomes 



rt.j=%(4V t , J+l + V, + l ,j+V,- l . } ). (5) 

In order to increase the accuracy of the relaxation 
method, we have expressed Vij in terms of the 
potentials of eight adjacent points. This is accom- 
plished by using Stirling differentiation formulas [2] 
for h\e points at equal intervals in both the r and z 
directions. The resulting formulas, including those 
special formulas required for points near the bound- 
aries, are listed in table 1. 

The basic nine-point formula used to relax the 
network is II, and its specialization to axial points 
is 12. Formula 13 was obtained from II assuming 
symmetry about the axis. Formulas 14 and 16 which 
were used for points near the left boundary are just 
special cases of the five-point formulas given in eqs 
(3) and (5). Formula 15 was obtained from II- and the 
relation 



- Figures in brackets indicate the literature references at the end of this paper. 
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Table 1. 


Formulas used for relaxation — Continued 




Number 


Potential 


Location 


Formula 
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Table 1. Formulas used for relaxation— Continued 



Number 


Potential 


Location 


Formula 
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Formula 17 was obtained from 12 and eq (6), and 18 
was obtained from 13 and eq (6). Formula 19 was de- 
rived by substituting into Laplace's equation the 
three-point formula for the derivative with respect 
to z and the five-point formulas for the derivatives 
with respect to r. Finally, Formula 110 was derived 
by substituting into Laplace's equation the three- 
point formulas for the derivatives with respect to r 
and the five-point formula for the derivative with 
respect to z. 

2.3. Over-relaxation Procedure 

The relaxation method is based on successive itera- 
tions of formulas II through 110 until the values 
assumed by the Vtj in two successive iterations 
differ by less than a preselected amount. Since the 
convergence is very slow, overrelaxation is used. 
The potential at point (i, j) after the mth iteration 
is then given by 



y.m.= ym-i+ w 

I, J I, J 



(v^-\ 



(7) 



where F* is the potential calculated from formulas 
II through 110. The quantity to is called the over- 
relaxation or acceleration factor, and its optimum 
value is problem dependent. 

To determine cu we use the method of Carre [3]. 
Twelve iterations are first performed using a value 
co =1.375 which is certainly smaller than the optimal 
value coo, and a quantity W is calculated from the last 
three iterations: 



W = 



y I y (m + i)_ y ( 

ikx \ y n V n 



m) 



i \Vn im) -v n (™-v 



(8) 



The next estimate of the acceleration factor is then 
calculated from 



(or 



l+\/l- 



JFio 2 



(9) 



Carre has shown the importance of not exceeding 
the optimum value a>o, and finds empirically that a 
value 



: <°'~7 (2 — &>z) 



(10) 



is close to coo but never exceeds it. The next twelve 
iterations are then performed using the acceleration 
factor (x) m and new estimates of W, o>/, and com calcu- 
lated from eqs (8-10). The entire process is repeated 
twelve times or until two successive estimates o>/, 
co/+i from eq (9) satisfy the inequality 



0)/ + 



o>i 



2- 



< 0.05 



-0)1 



(11) 



Succeeding iterations are then made with the last 
estimate of too from eq (9) and are repeated until 



Vf 



■ ym T \ 



W 



\-W 



(12) 



for all points on the mesh. The quantity e is the desired 
precision for the potentials. 

2.4. Program Organization 

The radius of the lens was divided into 40 mesh 
points, with 320 mesh points along the axis for the 
half-lens. Hence, the program must relax 12,800 points. 
In addition, a 5 X 100 matrix is relaxed in the gap be- 
tween the tubes, with 15 points overlapping the main 
matrix (see fig. 2). The time required to relax the com- 
plete network to a precision of e=10~ 8 was about 



B(l,40) 




(320, 40) 



A(l,l) 



F(320,l) 



FIGURE 2. Layout of mesh points showing overlapping of the main 
matrix with the auxiliary matrix used in the lens gap. Cross- 
hatched region is area of overlap. 
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TABLE 2. Axial potentials V(z)for the tWO-tube electrostatic lens 
2 = at the center of the gap. V\ = 0, V>= 1. 
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Present 
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Grivet- 




results 


et al. 


Bernard 




et al. 




Bernard 





.500000 


.500000 


.500000 


.500000 


.500000 


.500000 


.500000 


.025 


.467115 


.467206 


.467237 


.466896 


.466896 




.467089 


.05 


.434532 


.434703 


.434753 


.434101 


.434101 




.434462 


.075 


.402538 


.402789 


.402816 


.401911 


.401911 




.402393 


.1 


.371399 


.371715 


.371679 


.370599 


.370599 


.370606 


.371138 


.125 


.341351 


.341722 


.341566 


.340406 


.340406 




.340925 


.15 


.312592 


.313005 


.312673 


.311533 


.311533 




.311952 


.175 


.285279 


.285722 


.285159 


.284138 


.284138 




.284379 


.2 


.259528 


.259991 


.259146 
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.258334 


.258327 


.258327 


.225 
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.235881 


.234720 
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.234195 
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.25 
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.111851 
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.45 
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.070897 
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.6 
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.044536 
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.044009 


.044009 


.044012 


.040542 
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.027690 


.021183 


.027403 


.027403 


.027405 


.024331 


.8 


.017178 


.017141 


.014590 


.017006 


.017006 


.017003 


.014504 


.9 


.010642 


.010571 


.008658 


.010535 


.010535 


.010538 


.00861 1 
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.006586 


.006544 
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.006520 


.006520 


.006520 
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.004074 
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.000033 


.000015 
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.000013 
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.000005 
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.000008 


.000008 


.000003 


.000008 
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.000005 


.000005 
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.000005 


.000005 
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16 min on an IBM 360/65 computer, using double 
precision arithmetic (~ 15 decimal digits). 

2.5. Results 

In table 2 the potentials we obtain along the axis 
are compared with potentials obtained with other meth- 
ods. Our potentials are in agreement with Verster's 
calculated potentials [4] to within 0.3 percent. The 
systematic difference between our results and those 
of Verster is probably due to a difference in gap size, 
which in the present calculations is 0.1 D while Verster 
uses a zero gap. 

Also given are potentials calculated with the Grivet- 



Bernard approximation [6] for gap widths of 0.1 D 
and zero. The agreement with our results and also with 
the results of Verster is good close to the gap but 
becomes progressively poorer away from the gap as 
could be expected with the approximations used by 
Grivet and Bernard. Note that the systematic differ- 
ence between our results and those of Verster is in the 
same direction and roughly the same size as the differ- 
ences between the Grivet-Bernard potentials for the 
same gap sizes. 

We have also included the results of a calculation 
for the zero-gap case carried out by El-Kareh [5]. 
The agreement with Verster's calculation is extremely 
good. It should be noted that in deriving lens properties, 
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El-Kareh uses the Gri vet-Bernard approximation for 
the axial potentials. 

The Grivet-Bernard approximation assumes a linear 
potential variation across the gap, while our method 
of continuing the relaxation process into the gap shows 
that the actual potential variation in the gap differs 
substantially from linear, as shown in table 3. To 
assess the effect of the nonlinear potential variation 
in the gap, we made special calculations assuming a 
linear variation in the gap. It was found that the change 
in axial potential and lens properties was less than 
0.1 percent. 



Making the change of variables 



Table 3. 



Linearly interpolated potentials compared to values 
from relaxation in the lens gap 

, K. = 20. / = distance measured from edge of tube along inner diameter of the tube. 
.s = gap width = 0. 1 D. 



Us 


Interpolated 


Relaxation 





1.000 


1.000 


.125 


3.375 


3.987 


.25 


5.750 


6.315 


.375 


8.125 


8.441 


.5 


10.500 


10.500 


.625 


12.875 


12.559 


.75 


15.250 


14.685 


.875 


17.725 


17.013 


1.0 


20.000 


20.000 



Finally, very recent calculations of Read et al. 
[7] are given. For the zero-gap case, the same infinite 
series was evaluated by Read et al., as by Verster and 
by El-Kareh. The agreement is extremely good. For 
the finite-gap case, Read et al. use a new method in- 
volving a determination of charge distributions on the 
tubes which give agreement with the boundary 
conditions. The agreement with our results is within 
0.0005. 

3. The Ray Tracing Program 
3.1. Statement of the Problem 

Using the potentials obtained by relaxation we must 
solve numerically the Lorentz equation: 



F = -eE 



:i3) 



where F is the force on the electron, and E is the 
electric field which must be calculated from the 
potentials at the mesh points. We restrict to the non- 
relativistic case and to cylindrically symmetric elec- 
trostatic fields. In addition we consider only merid- 
ional rays, that is, rays lying in a plane through the 
axis of the lens. Equation (13) then reduces to the 
two equations: 



d^z_ 
dt 2 ' 

d 2 r_ 
dt 2 ' 



E z {r,z) 



■■ E r (r,z). 

m 



(14) 



2 2e ,2 

m 



dz 



v = 



d.T 

dr 
dr 



(15) 



eqs (15) become 



dz , 

du 1 , . 

-=--E z {r,z) 



(i6: 



dr 



■v{z, r) 



dv If/ -, 

-=-- Er (r,z). 



We express the initial conditions as 

z(to) =z 

(i)o = "° 

r(r ) =r 



(17) 



3.2. Method of Numerical Solution 

The numerical solution of the system of eqs (16) 
is of particular difficulty when a high precision is re- 
quired. After a number of attempts using previous 
methods [1] we have chosen to use the predictor- 
corrector method of Hamming [8] where error and sta- 
bility can be checked point by point. The important 
property of the classical predictor-corrector method 
is that it can be iterated until the required precision 
in the solution is reached. Hamming's method has the 
additional advantage of setting limits on the error which 
is accumulated along the integration. The method is 
stable in the sense that the difference between the 
numerical solution and the true solution decreases 
as the number of points of integration increases. 

The predictor-corrector method seems therefore 
ideal for the calculation of trajectories. It has, however, 
two disadvantages: (1) It is not "self-starting"; hence, 
several starting points must be calculated with a 
different method. (2) It is desirable that at least sev- 
eral points be calculated in each mesh. Since the points 
are calculated at equal time intervals h it may happen 
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that when the electron accelerates the stability criteria 
require the interval h to be reduced by a factor of 2. 
Then the predictor-corrector method needs, to proceed, 
points calculated at intervals h/2, half of which do not 
exist. The trajectory must therefore be reinitiated with 
a different method. 

Whenever required, the trajectory is reinitialized 
using the Runge-Kutta method. Reinitialization is 
only required when the electron moves toward the 
higher potential part of the lens. The effect of the 
occasional use of the Runge-Kutta method was checked 
by using the reversibility of electron trajectories. Sev- 
eral trajectories were calculated in both directions and 
found to coincide to within 0.01 percent. 

A negligible error is expected when the trajectory is 
begun, this process being always performed in a 
zero-field region. with 



n> = h • (t>t+2 P\ 
m-z = h - E z [ n -f- ~, Zi + -^ J 



ks= h • (iij + m 2 ) 
n :i = h ■ (vi + p-z) 
rriz = h - E z ( n + n-> , z\ + ki ) 
p 3 = h • E r ( r'j + n 2 , z-, + k-i ) 



3.3. Trajectory Initiation by the Runge-Kutta Method 

The numerical solution of eqs (16) uses the Runge- 
ivutta method to calculate the first seven points of 
the trajectory. Continuation of the solution then pro- 
ceeds with the predictor-corrector method. Whenever 
»he stability conditions fail to be satisfied at some point 
along the trajectory, the Runge-Kutta method is again 
applied to reinitialize the trajectory before continuing 
with the predictor-corrector method. 

The Runge-Kutta formulas utilized are as follows: 

zi+ 1 = zt + 7 ( Ao + 2A'i + 2A > + k-A ) 
o 



E Z = -EJ2 
E r =-Erl2 

E z and E, are obtained by differentiating with respect 
to r and z the Lagrange interpolating polynomial used 
to determine the potential at a given point from the 
values of the potentials at surrounding mesh points 
/(see sees. 2.5 and 2.6) . 

The time int erval h = At at the beginning is taken 
equal to l/(/? VF max ), thus assuring that the trajectory 
will be calculated for at least n points per mesh. Then 
h is reduced, if necessary, to satisfy the Hamming 
stability condition: 



where 



m+\ = Uj + - ( Wo + 2m i -h 2m > + m :{ ) 



i=n + 7 (n o +2rti+2rt2 + /*:0 




Vi+ 1 = vt + t ( po + 2pi + 2p> + p 3 ) 

ko = hiii 

/?<> = hvi 

m = hE z (ri, z%) 

p =hE r (ru zi) 



(18) 



ki = h - ( w/H- ^ m o) 



n t = h - lvi+j:po 



. no Ao 

E z \ri+ ~2' Zi ~2 

p x = h • E r [ n + y , zi + y 



fc> = A 



//,+ - mi 







ftA- < 0.4 






du 


+ 


dv 
dr 


+ 


du 


+ 





A= -^ + -^ 



(19) 



In this way, after the seventh point is calculated from 
the Runge-Kutta method, the predictor-corrector 
method is started with points spaced by an h which 
satisfies the stability condition. 

3.4. Trajectory Continuation by Hamming "Predictor- 
Corrector" Method 

The predictor-corrector formulas, as used in our 
program, are as follows: 

Predictor 



z\% = Zis + q h [2u, — Ui-t + 2ui- 2 ] 



*4°+V 



Ui-z-\-^h[2E z (n^i) 



E z (n-uzi-\) 
+ 2E z (zi- 2 Si- 2 )] 



rf?i = n_3 + 3 h [2vi - vt-i + 2t;«_ 2 ] 



^+i = *>*-3 + o A [2E Z ( nji ) - E z (n- u zi-\ ) 

+ 2^(r,_ 2 ,z,_ 2 )] 



(20) 
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Modifier 



A 2 6*,Mi 



*(o) =z (o) J.!!? / (on 

*i+l *i+l ' 221 

a}?, = < > i+Hf («i-«/ 0) ) 

112 



(21) 



(22) 



Corrector 



4iV , = 3 [9zi-z i -2 + 3h{u$ 1 + 2u i -u i -i)] 

O 



*4iV )= £ {9u,-u,_i + 3A[(£y)), + , + 2 (£,), 
r& ° = I [9r, - r t -2 + 3A ( !&>, + 2i>, - W _ , ) ] 

J i { + i l) =h9vi-v i - 2 + 3h[(Ey) i+1 + 2(Er)i 

-(£,),-,]} 



Ultimate Value [superscript / indicates the final 
values from eqs (22)] 
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h 2 b 2 -ik 2 k 2 



(l-hb-M) 



T<C, 

<c 2 



(24) 



where 



0-r 



Ar^la^-^-^j+i + lt/^-^-^li+i 



ti = 



A- 



*J = 



a£ 



+ 



+ 



a£V 



dz; 



"<+i w ;+i^ ^21 v '+ 1 < +1 



h+i r /+i^ 121 /+1 i+1 



(23) 



^t+i J+1 121 

The corrector formulae are iterated until the required 
precision in z, r, a, and t? is reached, that is, until the 
inequalities (24) are satisfied: 



du 

Cv=g(\z\ + \r\) 
C*=g(\u\ + \v\) 
g= maximum relative error 

Before calculating the ultimate values Zj+i, rj+i, 
^i + i, vt+u the interval h is checked against the sta- 
bility criterion (19). If satisfied, the next trajectory 
point is calculated with the same h. If not, h is divided 
by two until criterion (19) is satisfied. Starting with z,, 
//, ///, Vj seven points are then calculated with the 
Runge-Kutta method using the new h, and the solution 
continued. 

If the stability criterion (19) can be satisfied with an 
interval twice as large, h is doubled, and the next 
point with index i + 1 is calculated from the points of 
index i, i — 2, i — 4, and i — 6, instead of the usual /, 
i— 1, i — 2, and i — 3. 

3.5. Interpolation of Potentials 

Potentials between mesh points are calculated by 
Lagrange interpolation. Let n be the degree of inter- 
polation. Then 



f(x) = %C k (x)f(x k ) 



k=0 



(25) 



where 



C k (x) ■ 



(X — go) (X — X\) . . . {x — Xk-\){x — Xk+\) . . . (x—Xn) 
{x k — Xo)(Xk — Xi) . . . (Xk — Xk-l)(x k — Xk+l) . . . (x k — X n ) 



For the mesh we have used, the absciscae are uni- 
formly spaced with unit interval. Putting u — x — xo 
we have 



C k (x) = 



u(u-l) . . . (u-k+l)(u-k-l) . . . (u-n) 
(-l)"-*ifc!(n-*)! 
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(26) 



For our problem, the most reliable results were 
obtained with interpolation of order 4. 

The potential at point P(r, z) is calculated as follows: 
(1) The potential at each point indicated in figure 3 
with x and is interpolated from the potentials in the 
corresponding column or row. (2) Two values of the 
potential at P are then interpolated, one from the 
potentials at points labeled x and one from those 
labeled 0. The average; of these two values is taken as 
the potential at P. 

3.6. Calculation of Electric Field 

To find the two components of electric field, E z 
and E r , at any point, we differentiate the Lagrange 
interpolating polynomials [eqs (25) and (26)] with 
respect to z and r respectively. There results: 

C k (x) 

=2 u(u-l)...(u-i+l)(u-i-l)...(u-k+l)(u-k- \)...(u-n) 

£=0; if** 



(-l)"-*M(n-*)! 



(27) 



Here x represents either z or r as appropriate. Again 
rc = 4 was found to give the smoothest results and was 
used in the program. 

3.7. Accuracy of Results 

It is difficult to predict, a priori, the accuracies of 
the potentials and trajectories obtained from our 
program. The accuracy of the relaxation technique, 
using five-point formulas, is usually taken to be of 
order 1//V where N is the number of mesh points. 
Since we used the more accurate nine-point formulas 
we believe that the accuracy of the potentials is con- 
siderably better than 10~ 4 of the maximum potential. 

In the trajectory calculations the parameter g in eq 
(24) gives the minimum precision with which each 
point of the trajectory is calculated. Focal lengths 
obtained from paraxial trajectories calculated with 
g= 10~ 5 satisfy the relationship 



(28) 



to a precision of at least 10~ 5 , demonstrating that the 
estimate of precision is realistic. 

Even higher precision can be obtained by reducing 
g with a consequent increase in computer time. The 
use of g= 10 ~ 5 is sufficient to determine accurate 
first-order focal properties and third-order aberration 
coefficients for the two-tube lens. Results will be 
reported in separate papers [9], 



— X— X — 






p -*-* 



FIGURE 3. Points (x and 0) at which potentials arc interpolated to 
calculate the potential at an arbitrary point P(r,z). 
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